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Abstract 

Mixed-moment models, introduced in [141 129| for one space dimension, are a modification of the method of 
moments applied to a (linear) kinetic equation, by choosing mixtures of different partial moments. They 
are well-suited to handle such equations where collisions of particles are modelled with a Laplace-Beltrami 
operator. We generalize the concept of mixed moments to two dimension. The resulting hyperbolic sys¬ 
tem of equations has desirable properties, removing some drawbacks of the well-known Mi model. We 
furthermore provide a realizability theory for a first-order system of mixed moments by linking it to the 
corresponding quarter-moment theory. Additionally, we derive a type of Kershaw closures for mixed- and 
quarter-moment models, giving an efficient closure (compared to minimum-entropy models). The derived 
closures are investigated for different benchmark problems. 
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1. Introduction 

The full discretization of kinetic transport equations like the Fokker-Planck equation is in general very 
expensive since the discretized variable resides in X x S 1 2 x [0 ,tf] where X C R 3 and S 2 denotes the unit 
sphere in R 3 . Thus, the solution of the Fokker-Planck equation is very high-dimensional. 

A common approach to reduce the dimensionality is given by the method of moments mum. One chooses 
a set of angular basis functions, tests the Fokker-Planck equation with it and integrates over the angular 
variable, removing the angular dependence while getting a system of equations in x and t. Assuming now a 
specific form of the underlying distribution function, different approximate systems arise. Typical examples 
are the well-known spherical harmonics or Pjv models mm fl8] and their simplifications, the SPat [16] 
methods. These models are computationally inexpensive since they form an analytically closed system 
of hyperbolic differential equations. However, they suffer from severe drawbacks: The Pjv methods are 
generated by closing the balance equations with a distribution function which is a polynomial in the angular 
variable. This implies that this distribution function might be negative resulting in non-physical values 
like a negative particle density. Additionally, in many cases a very high number of moments is needed for a 
reasonable approximation of the transport solution. This is in particular true in beam cases, where the exact 
transport solution forms a Dirac delta. The entropy minimization M ,y-models m u in mug are expected 
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to overcome this problem since their closure functions are always positive. In many situations these models 
perform very well. Still they produce unphysical steady-state shocks due to a zero netflux problem. 

To improve this situation, half or partial moment models have been introduced in [10| I13j. These models work 
especially well in one space dimension, because they capture the potential discontinuity of the probability 
density in the angular variable which in ID is well-located. If however, a Laplace-Bcltrami operator is used 
instead of the standard integral scattering operator, i.e. scattering is extremely forward-peaked |261 . these 
half moment approximations fail [HIM]. 

To improve this situation a new model with mixed moments was proposed in HIM] which is able to avoid 
this problem. Instead of choosing full or half moments, a mixture of both is used. Contrary to a typical 
half moment approximation, the lowest order moment (density) is kept as a full moment while all higher 
moments are averaged over half-spaces. This ensures the continuity of the underlying distribution function. 

Since in one spatial dimension mixed moments perform very well, it seems reasonable to extend them to 
multi--D. It turns out, that here a mixture of full, half and quarter moments is necessary to derive the 
correct mixed-moment ansatz. 

Realizability is the fact that a vector of moments is physically relevant, i.e. that it is the moment of a non¬ 
negative distribution function. While in ID realizability theory for full moments [5j and mixed moments l2f)| 
is completely solved, it remains an open problem in higher dimensions. We will give a realizability theory 
for quarter and mixed moments of order 1 and derive, as in ID, a corresponding Kershaw closure which 
provides an analytically closed system of equations, in contrast to minimum-entropy models which requires 
the solution of a nonlinear system of equations. 

The paper is organized as follows: First, we will give a short introduction to the method of moments and its 
consequences for the Fokker-Planck equation. After setting up the basic notations in Section [2j we provide 
necessary and sufficient conditions of order 1 for realizability in the case of quarter moments and mixed 
moments in Section [3] Section [4] deals with the construction of Kershaw closures. Finally, we present some 
numerical tests for the derived models in Section [5] as well as conclusions and outlook in Section [6] 


2. Macroscopic Models 

We consider the Fokker-Planck equation 

d t i/j + fl ■ Vx^ + a a i/j = y + Q (2.1a) 

which describes the densities of particles with speed £ S 2 at position x £ X C R 3 and time t un¬ 
der the events of scattering (proportional to a s (f,x)), absorption (proportional to a a (t,x)) and emission 
(proportional to Q(f,x, fl)). The equation is supplemented with initial condition and Dirichlet boundary 
conditions: 

^>(0, x, Q) = i/’t=o( x , H) for x £ A', Q £ S 2 (2-lb) 

-0(f,x, Q) = x, fl) for t £ T, x £ dX, n • fl < 0 (2.1c) 


where n is the outward unit normal vector in x £ dX. 


Similar to m we assume that geometry, initial and boundary conditions are independent of the ^-direction, 
resulting in a solution ■)/> which is also z- independent. Therefore (2.1) can be reduced to x £ A' C M 2 . 


Parametrizing in spherical coordinates and taking the symmetry reduction into account we obtain 


12 = ^y/l — fj, 2 cos(ip), \Jl — n 2 sin(<p)^ =: (Aa 


Cly ) 1 


( 2 . 2 ) 
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where ip £ [0, 2tt] is the azimuthal and p € [— 1,1] the cosine of the polar angle. Then the Laplace-Beltrami 
operator on the unit sphere can be written as 


d 

dp 


dip' 


1 d 2 ip 


J.. I l 1 dfJ J ' l-^2 dip 2 


(2.3) 


Definition 2.1. The vector of functions b : S' 2 —> K n consisting of n basis functions bi, i = 0,.. .n — 1 of 
maximal order N (in ft) is called an angular basis. 

The so-called moments u = (u 0 ,... ,u n -\) T of a given distribution function ip are then defined by 

u = J hip dfl =: (hip) (2-4) 

s 2 


where the integration is performed component-wise. 


Equations for u can then be obtained by multiplying (2.1) with b and integration over S 2 : 

(b d t ip) + (bV x • flip) + (b cr a ip) = a s (b C (ip)) + (b Q) 


Collecting known terms, and interchanging integration and differentiation where possible, the moment sys¬ 
tem has the form 


d t u + d x {fl x hip) + d y (flyhip) + a a u = (b A n ip) + (b Q) 


(2.5) 


Depending on the choice of b the terms (fl x hip), (fl y hip) and in some cases even (bC (ip)) cannot be given 
explicitly in terms of u. Therefore an ansatz ip has to be made for ip closing the unknown terms. This is 
called the moment-closure problem. 

In this paper the ansatz density ip is reconstructed from the moments u by minimizing the entropy-functional 

U{ip) = (r 7 (V 0 > ( 2 - 6 ) 

under the moment constraints 

(hip) = u. (2.7) 


The kinetic entropy density ry : R —> R is strictly convex and twice continuously differentiable and the 
minimum is simply taken over all functions ip = ip(fl) such that Tl(ip) is well defined. The obtained ansatz 
ip = ip u , solving this constrained optimization problem, is given by 


V\i = argmin {(r](ip)) : (hip) = u} . 


( 2 . 8 ) 


This problem, which must be solved over the space-time mesh, is typically solved through its strictly convex 
finite-dimensional dual, 

a(u) := argmin (r]*(h T a.)) — u T a, (2.9) 

del" ' ' 


where ry* is the Legendre dual of ry. The first-order necessary conditions for the multipliers a(u) show that 
the solution to (2.8) has the form 


i’u = v, 


'* (b T «(u)) 


( 2 . 10 ) 


where ?y( is the derivative of ?y*. 
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This approach is called the minimum-entropy closure (21] . 

The kinetic entropy density rj can be chosen according to the physics being modelled. As in H3 nu, 
Maxwell-Boltzmann entropy 


vW = V’logW’) - ip 


( 2 . 11 ) 


is used, thus 77 * (p) = (p) = exp(p). This entropy is used for non-interacting particles as in an ideal gas or 

an ensemble of photons. 


A closed system of equations for u remains after substituting ip in (2.5) with ip u : 


d t u + d x (n x bip u ^ +d y (n.yhip u ^ + cr a u = y (bAaip u ^ + (b Q) 


( 2 . 12 ) 


Also note that using the entropy r](ijj) = \ip 2 the linear ansatz 


= b T a(u) 


(2.13) 


remains. If the angular basis is chosen as spherical harmonics of order N, (2.12) turns into the classical Pat 
model mmm- 


2.1. Angular bases 

This moment approach strongly depends on the choice of the ansatz ip and the angular basis b. In the 
following sections we will shortly derive the different angular bases for the models presented here. These 
bases will be generally collected in the basis-vector b(fl). If we need to further distinguish between the 
models, the corresponding symbols defined in the following sections will be used. If a result is independent 
on the choice of the basis we will just use b as symbol. 


2.1.1. Full moments 

The full-moment basis f at of order N consists of the tensorial powers of £2, i.e. (by abusing notation) 

fAr = (i, n, n ® n, n ® 3 ,..., n® N ) T . (2.14) 

The corresponding tensorial moment of order i is given by 

uW (2.15) 


which consists of the scalar moments 

1 ^, A) = i x ,i y >0,i x + iy = i. (2.16) 

Note that an equivalent system can be obtained by using the corresponding real-valued spherical harmonics 
of order N as angular basis 

2.1.2. Quarter moments 

We follow the approach in [12] where general partial as well as the special case of quarter moments in two 
space-dimensions are treated. The main idea is not to integrate over the whole sphere S 2 but over subsets 
of it. 
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Definition 2.2. For DCS 2 and ip £ L 2 (D,M.) we define its tensorial moment by 

u l D = (1 == <^V) D • 

As for full moments the corresponding components of the tensorial moments are given by 

= (D]fn l y y ip) D , i x ,i y >0,i x + i y = i. 

In this paper, D will be one of the following quarterspaces 


= 1 

'n | 

p £ [-1,1 ],<p£ , 

£++ — \ 

[n 

| p£ [-1,1], [0, |]}, 

S -H 

[n 1 

p £ [-l,l],<pe [tt, y]|, 

S + - = \ 

[ij 

1 p e [-1,1], G [^,2tt] j 

or halfspaces 







'n | 

p£[- i,i]^ e[ -|,|]}, 

Sf = \ 

[ n 

| M €[-i,i], v e[|,y]} 

1 

p£ [-1,1], ip£ [0, 7r] } , 

S- = {H I 

p £ [-1,1 },<p£ [tt, 2 tt]} . 

Note that for D = 


we have that u^j = 





(2.17) 

(2.18) 


For pure quarter moments, we will have D = Sij for i, j £ {+, — }. The corresponding basis for quadrant ij 
is then given by 


Ps, 


= l s • (1 


where • should be understood as multiplication with every component. Consequently the complete set of 
basis functions is p N = (ps ++ , Ps_+, Ps__, Ps + _) • 

2.1.3. Mixed moments 

As will be shown below the quarter-moment basis exhibits undesired properties when applied to the Laplace- 
Beltrami operator. This has inspired the works in mum where so-called mixed moments were developed. 
The main problem of the quarter-moment basis is that the ansatz ip is not continuous in f l which is necessary 


for the solution of (2.1) in one space-dimension. There, mixed moments where constructed by starting from 


a half-moment ansatz and demanding continuity of the ansatz (2.10) with respect to this basis. There, it 


suffices to choose a full zeroth-order moment and half-moments for all higher order moments [29] . 

The construction of mixed moments in two dimensions works in the same spirit. We start with the general 
quarter-moment basis p n and demand continuity of the ansatz ip u . 


Having (2.10) in mind, we obtain multipliers ay for every quadrant Sij. For example, for N = 1 we have 


V'u 


= v 


i (c 


(0,0) 


a 


111 



At the poles of the sphere (p = ±1) we have H = (0, 0) T , which implies that ip u is continuous only if 
agf = a(°’°) for all i,j £ {+, —}. Similarly it holds that along the quarter-space boundaries (i.e. ip = fcf, 
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k = 0,..., 3) some of the multipliers have to be the same (exactly those whose component of fl does not 
vanish on the boundary, e.g. a^ 1 ’^ = a^ 1 ) 0 ^ since at ip = 0 we have Vt x = 1 0). 

Accordingly we obtain the moments 

« (0 ’ 0) = (</>) = {n x n y if) D (2.i9) 

U s'P = (^x'0) 5 + Ug~ 0) = (^zV’)sy ( 2 - 20 ) 

4+° = ^y^s+ u s~ l) = {tov^Sy ( 2 - 21 ) 

for i = 1,..., N, i x + i y = i, i x ,i v > 0, and D € {<S ++ ,S_ + ,<S—,£+_}. 

In contrast to the one-dimensional setting one full moment, half moments for the basis functions contributing 
to either x or y direction, and quarter moments for the basis functions which contribute to both directions 
occur. 


Note that in the fully three-dimensional setting the decomposition in 3-direction has to be taken into account, 
leading to octants instead of quadrants. 

To embed this in the framework we choose our basis function m^r of order N as 


mjv = (l,f2 x l s +,...,f) x lg+jOels-,...,^ 1 5 


L V L Sy 5 ‘ • • ’ “y • 


L v L s~ > 


n x n y is ++ ,ftin y is ++ ,n x fi 2 v is ++) ..., 


n x n y i s+ _, fi 2 x n y i s+ _, n x n 2 y i s+ _) T 

=: (m S 2,m 5 +, m s -, m 5 +, m 5 -, m 5++ , m 5 , m 5 , m 5+ _) T . 


( 2 . 22 ) 


Formulas to efficiently calculate the appearing integrals in case of a linear ansatz can be found in Appendix 
|Appcndix A[ 


2.2. Moments of the Laplace-Beltrami operator 


All that remains to obtain a closed set of equations in (2.121 is to correctly evaluate (bAfj^ u ). This can 
be done using the formal self-adjointness of the Laplace-Beltrami operator, using 


(bA a ^) = (V’Anb). 


Due to this, the calculation of these integrals does not a priori depend on the choice of the ansatz but is 
true for every if. 


2.2.1. Full moments iN-basis 

Straight-forward calculations show that for i x ,i y > 0, i x + i y = i it holds that 

A. A',: <>;/ =-*(* + 1) + i x {i x - 1) h {iy - x ) 

Consequently, the corresponding moments are given by 

(fljf A n ip) = -i (i + 1) u ( bA H ) + ix ( ix _ !) u (C-2,i„) + iy _ x) 2)_ 


(2.23) 
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2.2.2. Quarter moments p n- basis 


We observe the following relations on the quadrant 5++ for i x ,i y > 0, i x + i y = i: 


Anls ++ = 


*'(¥>)+ <*'(! -<p) 


AqI 5. 


++^x 


1 — n 2 
rz —ix 2 

= V 1 -V 


2 " x +i x (i x - mi*- 2 i s++ - i x (i x + i)njf 15 . 


Anls ++ ^” = T 


'(f ~<p) + *v(*v~ W 2 


2 ** V(w- 


'5 


+ + 

— iy(iy + l)n“l 5+ 


An is,, = i s 


++ 

/ V <>v : + 


(-* (* + 1 ) + 4 (i x - 1 ) ni- 2 n]y +»„(*„- 1 ) n*-<>;/ 2 ) 


¥> - 2 ) ~ 5 (tO 


(2.24) 

(2.25) 


1 — /x 2 

Similarly, the corresponding quantities in the other quarter-spaces can be obtained. 

The moments of AnV’ include the evaluation of the microscopic values d^ip and ip at the quarter-sphere 
boundaries. Note that, similar to the one-dimensional case, the mass-conservation property of the Laplace- 


Beltrami operator (i.e. (A nip) = 0) is usually violated. This can be easily seen by summing up (2.24) over 
all quadrants, observing that d v ip at the angles k = 0,..., 3 remains. Also note that this quantity is 
not rigorously defined for minimum-entropy closures due to the discontinuity in the ansatz ip. 


2.2.3. Mixed moments m n- basis 
We observe for i x ,i y > 0 the following: 

An 1 = 0 

Afilofljf = — i x (ix + 1 ) 1 dHj + ix (ix — 1 ) 1 d ^ 1 x ~ 2 


±1 j)+a(y-ir) 

v'w* 

= iy (j"y + 1 ) 1 D^y H - iy (iy 1 ) l.D^ 
± H v =l 


y — vy \by ~r - l ; -l D* L y ■ °y K^y -*■; 

6(<p) + 6((p - it) 

" =1 


for D G {S+,S X } 

for D G {£+,£"}. 


The calculations for AqI s 3.V x r il’y are equivalent to those of the quarter-moment basis. 

Note that all these calculations are closure-independent. We therefore need to calculate u^± \ ul 0 ’ 1 *, 

S x c>y 

anc ^ ^he sem i" m i crosco Pi c quantities ip (p ,, fc|) dp for k = 0,... 3. 

3. Realizability 

In this section we will define precisely the concept of realizability. Furthermore we give necessary and 
sufficient conditions for first order quarter and mixed moment models. 

Definition 3.1 (Realizability). A moment vector u is said to be realizable with respect to basis b if there 
exists a non-negative distribution ip(Q) > 0 such that u = (b ip). The realizable set is defined as 

7Z b = {u € 1" | 3ip > 0 s.t. u = (b^)} (3.1) 


ip is then called a realizing distribution. 
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Note that (2.9) is solvable if and only if u € lZ h . 


A standard example for realizability are the realizability conditions of first order in 
Example 3.2. Let b = fi = (1, fl) T and u = (t/ 0,0 ), ul 1 1)^. Then u £ lZ fi if and 


iib 




the full-moment setting. 
only if m 

(3.2) 


3.1. Quarter moments 


In this section first order necessary and sufficient conditions for realizability of a quarter-moment vector will 
be given. This will be used later to derive the corresponding conditions for mixed moments. 

Lemma 3.3. For a vector of moments u = £ R 3 it is necessary and sufficient for the 

existence of a non-negative measure if which realizes u with respect to p$^ that 


and the normalized first moment 


satisfies <f> 


|i| 

Sij 


£ Sij. 


4 ' 


< «r 


== 


u 


|i| 


,( 0 . 0 ) 


(3.3) 


(3.4) 


Proof. We only prove the statement for <S ++ . The proof for the other quadrants works similarly. 
Assume that if) > 0 in 5++. Since ||fl|| 2 < 1 we obtain 


4 b 


W) 


s++ 


< (ll^ll 2 V’>< 


< = 4++ 


showing the necessity of (3.3). Since S2 £ 5++ we obtain Ll x , Q y > 0 implying that 4++ > ’ 4++ — 0- Together 
with (3.3) we have 4j + ^ <5++- 


To show the sufficiency of (3.3) we give a realizing distribution function. A possible (but not necessarily 
unique) candidate is given by 


ips ++ = 4++ )<5 ( n - <4++) 


(3.5) 


where (5 denotes the multi-dimensional Dirac-delta distributio 
is supported in <S++. Thus 


■fl 


if A'l + 


£ <S ++ the distribution function 


(if s++ ) = 4°;° + ) and ( n ^++) = 4°:?4 1 ! 


u 


lb 

■S++- 


Therefore, ips ++ is a realizing distribution for u under the given assumptions. 


□ 


'We assume for notational simplicity that 8 has mass 1 even on the boundary of integration. 
























3.2. Mixed moments 


With the knowledge of Lemma [373] we are able to provide the realizability conditions for mixed moments of 
order 1. 


Theorem 3.4 (First order necessary and sufficient conditions). For a vector of moments 

O x £>y £)y 

it is necessary and sufficient for the existence of a non-negative measure if which realizes u with respect to 
mi that 




(3.6) 


and Ugf\v,g+\ > 0. 

o~. o x 


Proof. We start again with the necessity of (3.6): 


Note that the vector of component-wise absolute values |fi| c := (|fl x |, \Cl y \) satisfies |||fi| c || 2 = ||fi|| 2 < 1. 
Therefore 1 — v T |fi| > 0 for every unitvector u £ S 1 , implying 


0< ((l-i' T |n| c ) if). 


This can be reformulated to 


v u := v 


T S. 


(1,0) (1,0) N 

u„ + — u 


1 * 


< u 


(0,0) 


(3.7) 


The left hand side will be extremal if v is collinear to u, i.e. v = yp=jp- Thus (3.7) implies (3.6). The 
sign-constraints on the moments follow again from the signs of f l x and Q y in the corresponding halfspaces. 


For sufficiency we give again a reproducing distribution. In the following we will make use that under (3.6) 
the reduced moment vector u is located in if 0,0 ^S ++ . Similarly the quantities 


n.o) (i,o) (o,i) _ (o,i)' 

Ji| | ■ u s+ u s~ s+ s~ 
<t>sl ■■= h ..(o,o) »•?- 


,( 0 , 0 ) 


hj £ { + > ~} 

will be in S, j . Reformulating (3.6) in terms of the normalized moments collected in 0^ we see that 


(3.8) 


V («>'■') := y(^»-^») 2 + (^'>-^« 


< 1 


(3.9) 


Also note that J\f ^0^^ = 05 - • We now embed the moments 0^. £ R 2 into the normalized mixed- 


moment space (which is a subset of R 4 ) in the following way. Define 

« + :=(C , -C-»’<f 1) -C , -°) T 

til :=(0,4">-4"',^«-4V>,o) T 

A'L := (o,*£°> - - 4V>) T 


(3.10a) 

(3.10b) 

(3.10c) 

(3.10d) 
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Figure 1: Linear interpolation between two realizability boundaries in the projected three-space 
(f>^+' > ) along an isoline of AM^ 1 '). The realizable set with respect to the quadrant 5_ + 
is plotted in grey. 


which then fulfils = ... = A/” 'j =■ c. Furthermore, any convex combination of those 

modified moment vectors along neighbouring quadrants (e.g. S ++ and S _|_) lies on an isoline of A/”(•), i.e. 

for ( G [0,1] we have 

A r (c4+ + + Cl - o4U) = c = Af + (1 - o4L) 

and analogously for the other half-space combinations. This is visualized in Figure [T] It can be shown that 
< fi ^ (and consequently u) can be written as a convex combination of the moments to . Indeed, 

defining 


C. ; = • „,I) e[MI srv 


4 V> 

Oy 


,(1^0) 
‘S X 


€ [ 0 , 1 ] 


(3.11) 


we see that 


<1 - + (> - wUl. = Kf - 

- (to'L + (i - <i)4'L = Kf - 


(o.i) 

Sv 


and finally 


0 |1| =C20' 1 | + (i-C 2 )0' 1 i. 


Since £ S^ we see that t/ 0,0 ) ^1, </> s can be realized (with respect to mi) by a distribution function 
with support in Sij, namely the quarter-moment distribution ipSij (realizing u*- 0,0 -* ^ 1, ^ with respect 


to Ps^ ) as given in equation (3.5) in Lemma 


3.3 


Therefore, due to the linearity of the problem, a non-negative realizing distribution for u is given by 

= (2 (CiVts+ 4 . + (1 - Ci)^s+_) + (1 - C 2 ) (Cl 4>S-+ + (1 - COV’S—) ■ (3- 12 ) 

□ 
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4. Kershaw closures 


A typical drawback of the minimum-entropy models defined by (2.101 is that the dual problem (2.9) cannot 
be solved analytically. The numerical solution, which has to be calculated at least once in every space-time 
cell, is challenging and expensive GO Dll- 


On the other hand, standard Pat models may give physically irrelevant solutions since they do not ensure 
positivity of the underlying distribution function. 


Due to this, Kershaw closures became recently a topic of increasing interest. They are constructed in such a 
way that they are automatically generated by a nonnegative distribution function ip. Therefore, the moment 
vector including the unknown highest moment is also realizable with respect to the basis of one order higher. 
Furthermore the flux function (ilhip^ is chosen to be exact (i.e. (flbip^ = (STb-0) where ip realizes the 

moment vector including the unknown highest moment) if u = Ui so = (b) is the isotropic moment. 

The last condition is also called the isotropic interpolation condition. 


In one spatial dimension for a full-moment basis, the Kershaw closure is also exact on the realizability 
boundary, because there the realizing measure on the realizability boundary is unique. This property is no 
longer true in higher dimension (see e.g. t 23J ) or for other models. However, it gives an idea of how to 
construct such a closure. 


The name of the closure is dedicated to David Kershaw who first proposed such an idea in [T9] . For an 
introduction into Kershaw closures in one space dimension we refer to [22] ■ The construction of Kershaw 
closures of first order as done in |29J requires realizability information of second order. With this it is 
possible to linearly interpolate between different parts of the higher-order realizability boundaries (choosing 
the interpolation parameter in such a way that the isotropic moment is interpolated). The resulting model 
is then cheap to evaluate because it is analytically closed (in contrast to minimum-entropy models). 

Unfortunately, we were not able to provide a closed second-order realizability theory for mixed or quarter 
moments yet which implies that we can’t identify the correct parts of this second-order realizability boundary 
for interpolation. However, it turns out that under some assumptions on <j> s . the second-order realizability 
information from the half moments in one space dimension are sufficient to define the Kershaw closure for 
quarter moments in a similar fashion. 

For mixed moments we choose a different approach to build the unknown second moment. Abusing the 
constructive procedure in Theorem |3.4| we are able to provide a closure for mixed moments by combining 
the quarter-moment closures accordingly. 


4-1. Quarter moments 


It was shown in [22] that, by assuming that the distribution function is symmetric around a preferred 
direction, the second moment of the Mi closure can be decomposed into 


J 2 I - 
<P S 2 ~ 


1-X 


■1 + 


(3x — 1) 052 052 


J 1 ! 

Vs 2 


(4.1) 


where X = X ( </>gi ) is the so-called Eddington factor. This implies that cp^i and (p^i are eigenvectors 


of . 


This is not necessarily true for all <pg but has to be true on a subset of the quarterspace. 
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Lemma 4.1. Let (fi s .. fulfil either 


( a ) 44 = °> (h) C’ 1; = °- ( c ) I fisf = C ] | or W ||^|| 2 = L 

TTien /or any < 4 ^ which is generated by a realizing distribution for we have that is an eigenvector 

of <j)g ]. with eigenvalue A satisfying 11</4 ,11 < A< Ill'll . 


Proof. We only prove the case «S»j = <S ++ . The other quadrants follow similarly. Assume that 0^ is 
generated by ip >0. Since fl x ,Cl y > 0 we have the following. 

(a) If (ps ’°^ = 0 it follows immediately that ip has support only on g> = | and n € [—1,1]. On this support 
we have f l x = 0 , which implies that 


<44 = 


S++ l 0 


0 0 


( 0 . 2 ) and 


/ 12 | ,| 1 | 4 0 \ / 0 \ ,( 0 , 2 ) ,\ 1 \ 

o 4?.v Ic/I s++ 

Similar to the one-dimensional case we observe that (keeping the support of -0 in mind) 

1 1 

4++ = (44)s ++ =2lr j ( x _ h 2 ) V’ d h < 27r J \A - dfj, = 

-1 -1 

and using Cauchy-Schwarz inequality 

l l 

“S++4+? = (44>s ++ Ws ++ = 47f2 J (! - 4) dfj, J ip dfi 

-1 -1 

> 471-2 f / \/l - d/j\ = 4°/ 1 + ) • 


Therefore the eigenvalue A := /4°’ 2 ' ) to the eigenvector </4 has to satisfy 


4 1 2 = 4°’ 1 ) 2 <a<4°’ 1) = 4 1 . 

2 *++ — — *++ 2 


(b) Analogously to (a) 


(c) If 4>Sif = 4°4 It follows that ip has support only on ip = j and n € [—1,1]. On this support we 
have Q x = fl y , which implies that 

/A(o,2) 

44 = (^fctf) Jfc® J , fading to 

41 41 _(4 2 4 4 2 ;°4 

4 + 4 + - ^o.o) 4,0) J (44 J ~ 2 ^ ++ 4 ++ - 
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(a) QM 1 eigenvector 


(b) Eigenvector deviation 


Ch + 



-Q- 


1 

0.8 

0.6 

0.4 

0.2 

0 





Figure 2: One eigenvector of the QM-l second moment and its minimal deviation from </>^| , measured 

in degrees. 


In this case we have 


u 


( 2 . 0 ) 

s++ 


(^)s ++ 




( 1 , 0 ) 

s++ 


while the lower bound is as before. Therefore the eigenvalue A := 2 (f>^’ 2 ' > to the 

eigenvector <j>g j has to satisfy 



= 2 u 


(i.o) 

5 ++ 


2 


< A < y/2’, 


t& 0) = 



(d) This follows exactly in the same way as for full moments, see e.g. [23]. Note that here, 


41 


= A = 


41 


= i. 


□ 


Numerical tests suggest that for the QMj model the first normalized moment is indeed always close (but 
not equal) to an eigenvector of the second normalized moment, i.e. remains small (where vi is 

the eigenvector of uj^. with smaller enclosing angle between itself and uj^.), see Figure jij This has been 
checked numerically using the QM-ptabulation strategy given in m 


Since the structure of eigenvectors of the QM : model is not immediately obvious to us, it seems reasonable 
to assume the simple form (4.1) for QK^ as well: 


41 = - 


4141 

a 2 --- - 2 

41 


(4.3) 
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where the coefficients aq(u),a! 2 (u) £ R>o have to be chosen accordingly. The eigenvalue associated to 
is given by A = oq + oq- 


41 


A short calculation shows that under this assumptions the bounds 

fulfilled for all (fig £ Sij, and not only on the subsets of Sij shown in Lemma 
holds that 


< A < 


4.1 


41 


have to be 


For every v £ Sij it 


4 v = ( n ) < (ttip) Sij = u s. 


|i| 


<i 


where the last inequality is meant component-wise. This implies for u = 




K 1 !, 


that 


u|o u = Xv < ul 11 , 

tJi-j On' ' 


and consequently A < 



2 


1 2 1 

The lower bound on A follows in the same way as for full moments observing that 
be positive semi-definite. 




T 


lias to 


The construction of the Kershaw closure in one spatial dimension linearly interpolates between the upper 
and lower boundary values for the second moment in such a way that the isotropic point is correctly hit. 
We define this interpolation not on the second moment but on its eigenvalues. The isotropic moment 
vector is given by 4’°+ = 4+^ = \ (similarly for the other quadrants with appropriately changed sign), 

,(2,o) _ ,(o,2) _ i _, + ,(i,i) 

^5++ — PS + + — 3 anCl ^5 ++ — 37T • 

We therefore define 


A — Oq + 02 — C 


41 


+ (i-C) 


41 


where £ = ( (A|sj + ) G [0,1] has to be chosen accordingly. Plugging in the isotropic point we get that 

c(|.|) = - 2 ,V(^ “ 0 - 7801 - 

Although other choices are possible we furthermore assume that oq + a-i is rotationally symmetric, which 


implies that £ = £ ^ (fi'g,, J. Unfortunately, the simplest choice £ = £ (|, |) leads to a moment system 
which is no longer hyperbolic. In fact, some of the eigenvalues of the flux Jacobian are imaginary, especially 
in those regions where the QM-l eigenvectors deviate strongly from (fig (compare Figure 2). We fixed this 

problem by choosing £ to be linear in 




i.e. 






The proof of Lemma 4.1 reveals more information about the choice of aq. In situations (a), (b) and (d) (i.e. 
on the realizability boundary), oq has to be zero, while it has to be 1 — ^ at the isotropic point (recall 
that there ^ due to the off-diagonals of the second moment). The simplest function fulfilling these 

two conditions is given by 
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Figure 3: Components of the QK : second moment on <S ++ . 


a x 







The corresponding closure relation is shown in Figure [3] Despite the fact that the assumption on the 
eigenvectors for QM-l is not true for all (j)g . the resulting second moment for QK^ is very close to the 
minimum-entropy one. This is shown in Figure [4] 

Since we are not able to provide a closed second-order realizability theory we cannot conclude immediately 
that this closure is in fact realizable with respect to P 2 . We therefore have to go a different way and, as 
before, provide a realizing distribution for this special case. We only prove the case Sij = £> ++ and drop the 
lower index <S ++ for readability issues. We make the ansat^] 

i> = u^' 0) ^2c L ±S(n - ± p t ±v t )), (4.4) 

where l £ {1,2} and the parameters c L ± 7 p L ± £ R and v t £ R 2 still have to be determined. Without loss of 
generality we assume u^ 0,0 ^ = 1 for brevity. The moments of this ansatz are 

1 = Ws++ =X^±’ 

0 111 = (D^) 5++ = + X^ Cl + Pt + “ 

0 |21 = (nn T ip) s++ = 0 | 1 | 0 |11 +^{cl+Pl+ -c L -p L -)\ L 4> w 

+ ( Cl +^?+ + D-P?-) v t vf. 


2 We acknowledge G. Alldredge for presenting us this ansatz. 
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Deviation in 


Deviation in ^j 1 ’^ 



Figure 4: 


Deviation of 


<?5 ++ 


and 


A2,0) 

0S++ 


for QM-l and QK 1 . respectively. 


Hierarchically inserting the equations into the one of higher order we see that the coefficients c L ± , p L ± fulfil 
the system of equations 


1 = H Ct± 

0 — Cl-\-Pl-\- C-l—Pl— 

"V C L+Pl -f 4“ C L—Pi — 


tG {1,2}, 
4 6 {1,2}. 


and (A t ,v t ) is the t-th eigenvalue-eigenvector pair of 0^ — 0^0^ , i.e. 

i _L 


Vl = 


V 2 = 


0 

0 


Ill- 


Ill 


0 


111 


0 


111 


Ai — ai 


A 2 = oti + — 


0 


111 


(4.5a) 

(4.5b) 

(4.5c) 


Since ip has to be supported in 6> ++ and should be nonnegative we also get the inequality-constraints 

c L ± > 0 (4.6a) 

0 |1! ± p t ±v t G £++ (4.6b) 


Combining (4.6a) and (4.5b) shows that sign (p,,+) = sign (p t _), thus, without loss of generality, the quarter- 
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1 



111 T 

Figure 5: Visualization of the coefficients p L ± for the example = (0.4,0.6) . The length of the arrows 
represent the maximal value for the respective coefficient. 


sphere constraints (4.6b) can be rewritten as 


0 < pi+ < min 


0 < pi- < min 


U'> 

Us? 




0 < P2+ < 1 - 


0 


111 


— : P 2+ 


0<p 2 -< 


0 


111 


= : P2- 



= ■ Pl + 



=' Pl- 


A visualization of these bounds can be found in Figure [5] for a specific example. The min function ensures 
that the first intersection with either the coordinate axes (pi+ in the example) or the norm-bound (pi~ in 
the example) is taken. For p 2 ± there is no distinction necessary since we can only intersect with one of 
them. 


Solving (4.5b) and (4.5c) for c L ± we get 


c t ± = 


V 

Pl± (p t + + p«,-)' 


Plugging this into 


(4.5a) we obtain an equation for p L ±\ 


Mpi+pi- + A 1 P 2 +P 2 - — P 1 +P 1 -P 2 +P 2 - 


(4.7) 


We found that 


p i_ = mm 


Ao, 1 ) 

p s ++ 

Abo) 

p s ++ 


0 


111 


f - 


0 


111 


P 2- = 


111 


0 

2A,. r , 

Pl+ — -j 1 € {1, 2} 

Pi- 


solves (4.7) while fulfilling (4.6b). 
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We want to prove this exemplarily for P 2 +- After some symbolic simplifications we see that 


P 2 + = 


2 y/2 (4 - 7T) (y/2 + 1) 


37r 


4> 


lb 


0 


m 


> o. 


With this it immediately follows that 1 — 




id 


- P2+ > 0. 


Note that the above choice of the coefficients p,± is not unique. Since the corresponding values for c L ± 
are also non-negative, (4.4) is a non-negative distribution function with support in <S ++ which realizes the 


desired moments. Therefore, the QK^ closure is realizable. 


4-2. Mixed moments 

Since the mixed-moment setting is even more complicated than the quarter-moment one, we want to follow a 
different approach to construct an approximate closure. As has been shown in Theorem |3.4| it suffices to have 
a non-negative quarter-moment distribution ips, ;l to construct a mixed-moment reproducing distribution. If 
these distributions fulfil the isotropic moment interpolation condition for the quarter-moment basis, by 


linearity, the resulting mixed-moment distribution (3.12) will fulfil it as well for the mixed-moment basis. 


Such a distribution function for quarter moments has been found in (4.4). Thus setting </>£ . as in (3.8) will 


provide the desired mixed-moment distribution function if mi using the interpolation defined by (3.11). 


A closure can now be generated by calculating the second moments with respect to those quarter-space 
distributions which have support in the corresponding set D. For example, 

m (2 ; 0) = (f+ (n 2 x i/; mi ) = C2Ci«s+? + C 1 - 


where is defined as in (4.3) with substituted by cj>g 


'S+- 
|i| 


Remark 4.2. Numerical experiments suggest that this closure is also (weakly) hyperbolic. No imaginary 
eigenvalues of the flux Jacobian have been found on a finely discretized grid of the realizable set. However, 
at the isotropic point both flux Jacobians have three times the eigenvalue 0. Still, the transformation matrix 
has full rank, i.e. the geometric multiplicity of the eigenvalue 0 is 3. 

Remark 4.3. Any other choice of'ifsij which fulfils the desired properties with respect to the quarter-moment 
basis will be a feasible choice as well. Exemplarily, the quarter-moment minimum-entropy closure QM 1 shall 
be mentioned. An efficient implementation using tabulation is given in 


Note that this closure is not an approximation of the MMj model and might behave differently. This is 


shown exemplarily in Figure 


where the second moment for the MMj model and the Kershaw closure 


is plotted against and with = — 4>^+' 1 = 1 fixed. The corresponding shape of the Kershaw 

Ox Oj Oy iOy 

closure is very different to those of the MMj model. 


4-3. Treatment of the Laplace-Beltrami operator 

A clear drawback of the mixed-moment Kershaw closure is that the moments of the Laplace-Beltrami op¬ 


erator (see Section 2.2.3) require the evaluation of the ansatz function. When using the QK]^ distribution 
to realize the mixed-moment closure, the ansatz is a linear combination of Dirac deltas, which cannot be 
evaluated in a feasible way. 
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(a) Kershaw (b) MM] 

MM 4 °'m- 0,155 °'| 16 - °-( 6 5 0.17 0.175 



One could follow the authors in [14] and replace the ansatz in the Laplace-Beltrami moments by the corre¬ 
sponding polynomial distribution (2.13). Then, we have (miA^i/)) ss S u with 


S = 


( 0 
do 
—do 
do 
\ —do 


0 

do 

d 2 

do 


0 

d 2 

do 

-do 

do 


0 

do 

—do 

do 

d 2 


0 \ 
-do 
do 

d 2 

do 


do — — 
do = 


- 1 


7T — 4 
3 7T (-7T — 3) 


0?2 — — 
do = 


(2 7T — 4) ( 7 r — 4) 
3 7T 


(2 7T - 4) O - 4) 2 

3 7T 1 

(2 7r — 4) (tt-4) ~~ 2' 


(4.8) 


Unfortunately, this closure may lead to an undesired behaviour of the moment-system (2.12), since its so¬ 


lution may leave the realizable set (this can be easily shown even in one spatial dimension with the ansatz 
given in El)- This follows from the fact that the tangent field spanned by 5u may point outside of the 
realizable set at those parts of the boundary where the polynomial ansatz (2.13) is negative. 


Consider for example u = (1,1,0,0, 0) . Then we have that 


S'u « (0, -2.813,0.813,1.813, -1.813) T , 
pointing outside the realizable set in the third component (recall u^' 1 < 0). 


We therefore modify this approach slightly and replace the Dirac ansatz with the tabulated QM]^ distribution. 
Since the fluxes for QM-] and QK 1 are very similar, it can be expected that the error done with this 


replacement is negligible. Furthermore, since the ansatz is positive, the moment-system (2.12) will behave 
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as expected. Indeed, calculating the approximation for the previous example we get approximately 

(m 1 A n i/'u) « (0,-2,0 ,d,-d) T , d -> oo, 

which points inside of the realizable set. Note the fourth and fifth component which converges to ±oo. This 
is consistent with the original problem since on the realizability boundary, the exponential ansatz converges 
to a combination of Dirac deltas. 


5. Numerical results 

We present the derived models in several benchmark test-cases. The numerical discretizations are obtained 
with a two-dimensional generalization of the high-order realizability-preserving discontinuous-Galerkin scheme 
given in [2] for the Kershaw models and a generalization of the realizability-preserving kinetic scheme given 
in [301 for minimum-entropy models. Those generalizations are in principle straight-forward but might 
be topic of a follow-up paper. The spatial and temporal order is four where roughly 10000 rectangular 
(discontinuous-Galerkin scheme) and 20000 triangular elements (kinetic scheme) in space were used. The 
reference solution is discretized with the second-order scheme given in mmm- 


5.1. Line Source 

The line-source test is a Green’s function problem, where a pulse of particles is emitted in an infinite medium 
[6]. It has been widely investigated for isotropic scattering in [15]. We choose the following set of parameters 
for this problem, smoothing the initial Dirac delta in space: 

• Domain: X = [-§, |] x [-±, §] 

• Final time tf = 0.45 

• Parameters: a a = Q = 0, cr s = 1 

exp (- x 2 +y 2 ) 

• Initial condition: ^i t=0 (x, Cl) = max(—^— 2 2 g2 ' , 10 -4 ) with a = 0.03 


Boundary conditions: Isotropic in Cl, consistent with initial conditions. 


We show several model solutions in Figures [7] and [8j We observe that convergence is achieved quickly with 
increasing moment order N for M^v models which can be observed as well in one spatial dimension. This is a 
huge contrast to the integral-scattering operator used in HE where only slow convergence can be observed. 
Note that, by construction, the MMi model has no rotational symmetry. The main directions of propagation 
follow the half-spaces, i.e. along the cartesian axes. Still, it performs significantly better than Mi (especially 
along the diagonal cut shown in Figure 5.8b), but worse than M 2 which almost converged to the reference 


solution. We observed that the naive implementation of higher-order mixed-moment closures like in the MM 2 
model are numerically instable due to errors in the non-linear closure of the Laplace-Beltrami operator (see 


Section 2.2.31. 


Figures 5.7e and 5.7f show the MKi model with the two different closures for the Laplace-Beltrami operator 


given in Section 4.3 The polynomial closure (4.8) has a slightly lower peak speed as the tabulated closure 


but more interesting is that the latter is much closer to the MMj solution (Figure [5.7d[ ). This shows, that in 
contrast to the one-dimensional situation in HU: the choice of the Laplace-Beltrami closure has a significant 
effect on the solution of the mixed-moment system. 


Since a s > 0, the QKj model cannot be applied here. 
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Figure 7: Local particle density z/ 0,0 ) of some models at t = 0.45 in the line-source test case. 
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(a) y = 0 


(b) x = y 


—■•— Reference —•— Mi —o— M 2 




(c) y = 0 (d) x = y 

—•— Reference —MMj —MKi poly MKi 




Figure 8: Local particle density id 0,0 ) of some models at t = 0.45 in the line-source test case, horizontal and 
diagonal cuts. 
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5.2. Two Beams 


This is the typical situation where the Mi model completely fails and is an often-used benchmark problem 
in fD (see e.g. [HU2S])- Instead of using opposing beams we investigate two orthogonal beams hitting each 
other in a void. We show two different variants of this. The first one is given by 


• Domain: X = [0,1] x [0,1] 


• Final time: tf = 1.2 


• Parameters: a s = a a = Q = 0 


• Initial condition: ^; t _ 0 (x, ft) 

• Boundary conditions: 


1CT 4 

47T 


V> b (t,X, ft) 


100 

47T 


exp 


exp 


(-=*£) 

^ M 2 + (v ~ f) 


10 -6 


if ar = 0, y e [0.45,0.55] 

if y = 0, X e [0.45,0.55] 
else 


where cr 2 = 0.05. 


Figure [9] shows the solutions of full- and mixed-moment minimum-entropy models up to order iV = 3, 

In the second version, the whole setup is rotated clockwise. The two beams sit in the upper-left and lower- 
left corner, respectively, both pointing to the center of the domain. Due to the rotational symmetry of the 
full-moment models and the reference solution, only the MKi closure is shown in Figure El It is visible 
that in contrast to the first version of this problem (compare Figure 5.9h I no shock is produced when the 
beams hit. 


6. Conclusions and future work 

We developed a two-dimensional extension of mixed moments, leading to a generalization of the original 
minimum-entropy MMjv model, proposed in [14[ l29j. Additionally, we provided a first-order realizability 
theory for mixed as well as quarter moments. Since the numerical inversion of the general minimum-entropy 
moment problem is very expensive we developed approximate, but still realizable, closures, the so-called 
quarter- and mixed-moment Kershaw closures. They have approximately the same cost as the corresponding 
polynomial closure while being realizable and non-linear, better adapting to the correct eigenvalues of the 
true Fokker-Planck solution. Although first-order mixed moments completely resolved the zero-netflux 
problem of the Mi model, a first order approximation in two dimensions is in general not enough, as has 
been shown in the two-beams test-case. Under some assumptions on the alignment of the beams this 
drawback can be completely removed. On the other hand, going to higher-order approximations, like Mjv 
or MMjv with TV > 2, a much better approximation can be done. Thus it seems reasonable to extend the 
concept of Kershaw closures for either K^v or MK^v for TV > 2 to obtain a cheap and accurate approximation 
of the Fokker-Planck equation. Furthermore, more research is necessary to obtain better-quality high-order 
numerical solutions of Kershaw closures. 
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Reference Mi M 2 



Figure 9: Local particle density u (°’°) of some models at t = 1.2 in the two-beams test case. Colorscale cut 
at u(°’°) = 3. 
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(a) y = 0.5 


(b) x = y 


—9— Reference — Mi —o— M 2 * M 3 



(c) y = 0.5 (d) x = y 

—9— Reference — MM] — 0 — MM 3 —9— MKj 




Figure 10: Local particle density of some models at t = 1.2 in the two-beams test case, horizontal and 
diagonal cuts. 
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t = 0.85 


t = 1.7 



x 
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Figure 11: Local particle density i/ 0,0 - 1 of the MKi model at different times in the rotated two-beams test 
case. 


Appendix A. Calculation of partial-moment integrals 


< n iw>s ++ = 

< n w> 5 _ + = 
<w> s __ = 
< n w> 5+ _ = 
( ni ) s + = 

where B ( z , w) is the beta function. 




^3; -f- 1 fy 1 1 


4^ B G > i + ^ (ix+i, ')) B 


2 ’ 2 

Za; “hi Zy “h 1 


(- 1 ) 




® ( o ’ 1 + O + iy) ) B 


2 ’ 2 
i 


1 iy + 1 


+ 1 *y + 1 


(-1 ) iy (\ 1 

2- 1+ 2( <x+< «') 

(- i ) < ( n x> 5 - = ( n i) s .t = (-i) < ( n i) s .T = 7 


2 ’ 2 
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